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We present numerical A^-body simulation studies of large-scale structure formation. The main 
purpose of these studies is to analyze the several models of dark matter and the role they played 
in the process of large-scale structure formation. We analyze in this work a flat cold dark matter 
dominated model known as the Santa Barbara cluster. We compare the results for this model 
using the standard Newtonian limit of general relativity with the corresponding results of using the 
Newtonian limit of scalar-tensor theories. An specific model is the one that considers that the scalar 
field is non-minimally coupled to the Ricci scalar in the Einstein-Hilbert Lagrangian. Comparisons 
of the models are done showing results of rotation curves, density profiles, and velocity dispersions 
for halos formed at z=0. We analyze, in particular, the Santa Barbara cluster and its possible 
equation of state. 

PACS numbers: 95.30.Sf; 95.35.+d; 98.65.-r; 98.65.Dx 



I. INTRODUCTION 

The Santa Barbara (SB) cluster model was introduced by Frenk et al. [1] in order to study in a systematic way a flat 
cold dark matter (FCDM) dominated universe using a variety of numerical codes. The main goal of this comparison 
was to asses the reliability of cosmological simulations of clusters in one of the simplest astrophysical relevant case. 
They compared the images and global properties of the cluster obtained by the different numerical codes. Heitmann 
et al. (2005) [2] analyze again the cluster with others new numerical codes and with a similar purpose and now it has 
become one of the standard cases of study. Given that this is a standard case to test we have decided to analyze it 
in the framework of the scalar-tensor theories (SST) of gravity. 

Scalar fields have been considered as one of the best possible ways to modify gravity. The work by Nordstrom, 
published before general relativity, formulated a conformally flat scalar theory of gravity Q , and finally, the scalar 
field role in gravity has been stablished since the pioneering work of Jordan, Brans, and Dicke[4, 5]. Nowadays 
they are considered as a mechanism for inflation [g.]; the dark matter component of galaxies [7]; the quientessence 
field to explain dark energy in the universe Q. The main goal of this work is to study the large scale structure 
formation where the usual approach is that the evolution of the initial primordial fluctuation energy density fields 
evolve following Newtonian mechanics in an expanding background [9]. The force between particles are the standard 
Newtonian gravitational force. Now, we will see that we can introduce the scalar fields by adding a term in this force. 
This force will turn out to be of Yukawa type with two parameters (a, A)[lo']. For so many years this kind of force, the 
so called fifth force, was thoroughly studied theoretically [11] and many experiments were done to constrain the Yukawa 
parameters [Ijl. We have also been studying, in the past years, the effects of this kind of force on some astrophysical 
phenomena[ifl, [H, 13, [H and in cosmological simulations |17|. The Yukawa force comes as a Newtonian limit of 
a scalar-tensor theory with the scalar field non-minimally coupled to gravitation [l^ although other alternatives can 
be found [19]. 

Our general purpose is to find the role these scalar fields play on the large scale structure formation processes. In 
particular, in this work we present some results about the role scalar fields play on cosmological simulations that form 
the SB cluster. 

We start by discussing a FCDM model and the general approach in A/'-body simulations (See Bertschinger[20] for 
details). Then, we present the modifications we need to do to consider the effects of a static scalar field and we 
show the results of this theory for the FCDM model that form the Santa Barabara cluster P, To perform the 
simulations we have modified a standard treecode the author has developed [2l| and the Gadget 1 [22[ (see also 
|http: //www. astro, inin. mx/mar) in order to take into account the contribution of the Yukawa potential. We finish 
this paper by discussing how we can obtain the equation of state for a dark matter halo in the framework of general 
relativity and its Newtonian limit. 
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II. EVOLUTION EQUATIONS FOR A CMD UNIVERSE 
A. General Scalar-tensor theory 



The Einstein equations for a typical scalar-tensor theory with a massive scalar field non- minimally coupled to the 
geometry are given by 



and the scalar field equation 



(t)V' - 2V 



[SttT- 



(1) 



(2) 



3 + 2cj 3 + 2cj 

where {)' = Here ^^^^ is the metric, R is the Ricci's scalar, i?^^^ the Ricci's tensor, T^^^ is the energy-momentum 
tensor, and u;{(j)) and V{(j)) are arbitrary functions of the scalar field (j). We will not consider a cosmological constant 
contribution in this work. 



B. Newtonian approximation of STT 

The study of large-scale formation in the universe is greatly simplified by the fact that a limiting approximation 
of general relativity, Newtonian mechanics, applies in a region small compared to the Hubble length cH~^ i^^o^ ^ 
3000/i~^ Mpc, where c is the speed of light, Hq = lOOh km/s/Mpc, is Hubble's constant at this epoch and h ^ 0.7), 
and large compared to the Schwarzschild radii of any collapsed objects. The rest of the universe affect the region only 
through a tidal field. The length scale cH^^ is of the order of the largest scales currently accessible in cosmological 
observations and H^^ ^ 10^^ yr characterizes the evolutionary time scale of the universe. 

Therefore we need to describe the STT theory in its Newtonian approximation, that is, where gravity and the scalar 
fields are weak (and time independent) and velocities of dark matter particles are non-relativistic. We expect to have 
small deviations of the metric with respect to Minkowski metric and of the scalar field around the background field, 
defined here as {(j)) and can be understood as the scalar field beyond all matter. If one defines the perturbations 
^ = (j) — {(j)) and hfj^y = g^y — rj^y^ where rj^y is the Minkowski metric, the Newtonian approximation gives [18| 

Rm = h^h^o = T^47rp - \v^ , (3) 

V^^ - mlp4) = -Snap , (4) 

we have set {(/)) = (1 + a)/GN and a = 1/(3 + 2uj). We are considering that the infiuence of dark matter is due 
to a boson field of mass msr governed by Eq. (|4]), that is the modified Helmholtz equation. Equations (|3]) and (|4]) 
represent the Newtonian limit of a STT with arbitrary potential V{(j)) and function u;{(j)) that where Taylor expanded 
around {(j)). The resulting equations are then distinguished by the constants On (the local gravitational constant), 
a, and A = hp/msFC. Here hp is Planck's constant. 

Note that Eq. (|3]) can be cast as a Poisson equation for i/j = (l/2)(/ioo + 4^/{^))^ 

VV = 47r-^p. (5) 

1 + a 

The next step is to find solutions for this new Newtonian potential given a density profile, that is, to find the 
so-called potential-density pairs. General solutions to Eqs. (|4]) a nd ([5j) can be found in terms of the corresponding 
Green functions, and the new Newtonian potential is (see ^0|, [iJ] for details) 



= -hoo = — / dr, 
2 1 + <^ J 



■ „.^/,, /('f-"-;-'^\ B.C. (6) 



The first term of Eq. (|6]), given by ?/^, is the contribution of the usual Newtonian gravitation (without scalar fields), 
while information about the scalar field is contained in the second term, that is, arising from the infiuence function 
determined by the modified Helmholtz Green function, where the coupling uj{a) enters as part of a source factor. 
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C. Cosmological evolution equations using a static STT 



To simulate cosmological systems, the expansion of the universe has to be taken into account. Also, to determine 
the nature of the cosmological model we need to determine the composition of the universe, i. e., we need to give the 
values of fli for each component z, taking into account in this way all forms of energy densities that exist at present. 

If a particular kind of energy density is described by an equation of state of the form p = wp, where p is the pressure 
and w is di constant, then the equation for energy conservation in an expanding background, d{pa^) = —pd{a^)^ can 
be integrated to give p oc a~^^^~^'^\ Then, the Friedmann equation for the expansion factor a{t) is written as 

i 

where k characterizes the geometry of the universe {k = for a flat universe), and Wi characterizes the equation of 
state of specie i. Here, = pi/ pc^ with pc = 3H'^ /SttGn- The most familiar forms of energy densities are those due 
to pressureless matter with Wi = (that is, nonrelativistic matter with rest-mass-energy density pc^ dominating over 
the kinetic-energy density pv'^/2) with ^dm ^ 0.22 and radiation with Wi = 1/3 and ^ 2 x 10~^. The density 
parameter contributed today by visible, nonrelativistic, baryonic matter in the universe is Qb ^ 0.04. There is also 
another main component, the energy density associated to a cosmological constant, ^0.74 with an equation of 
state such that Wi = —1. In this work we will consider a model with only one energy density contribution. One which 
is a pressureless and nonbaryonic dark matter with contribution given by Qdm that does not couple with radiation. 
The above equation for a{t) in this case becomes 

^^9^ /<^o\^ k 



p) -4 (8) 



Here, we employ a cosmological model with a static scalar field which is consistent with the Newtonian limit given 
by Eq. (|6]). Thus, the scale factor, a(t), is given by the following Friedman model (see a more general case in [17]), 



^DMO , ( ^ ^DMO 



1 + a \ 1 



(9) 



where ^dmo is the dark matter density evaluated at present, respectively. We notice that the source of the cosmic 
evolution is deviated by the term 1 -\- a when compared to the standard Friedman-Lemaitre model. Therefore, it is 
convenient to define a new dark matter density parameter by ^^^If = ^dm/(1 + Q^)- This new density parameter is 
such that always ^^^m ~ which implies a fiat dark matter dominated universe, and this shall be assumed in our 
following computations. For positive values of a, a fiat cosmological model demands to have a factor (1+a) more 
energetic content {^dm) than in standard FCDM cosmology. On the other hand, for negative values of a one needs 
a factor (1 + a) less Qdm to have a fiat universe. To be consistent with the CMB spectrum and structure formation 
numerical experiments, cosmological constraints must be applied on a in order for it to be within the range (—1,1) 
(23|, M I25|,l26i. 

In the Newtonian limit of STT of gravity, the Newtonian motion equation for a particle i is written as 

where x is the comovil coordinate, and the sum includes all periodic images of particle j, and Fsri^^ A) is 

FsF{r, a. A) = 1 + a (l + ^) e'^^^ (11) 

which, for small distances compared to A, is Fsri^ < A, a. A) 1 + a (l + J) and, for long distances, is Fsf{'^ > 
A, a. A) ~ 1, as in Newtonian physics. 

We now analyze the general effect that the constant a has on the dynamics. The role of a in our approach is as 
follows. On one hand, to construct a fiat model we have set the condition ^^2f ~ which implies having (1 + a) 
times the energetic content of the standard FCDM model. This essentially means that we have an increment by a 
factor of (1 + a) times the amount of matter, for positive values of a, or a reduction of the same factor for negative 
values of a. Increasing or reducing this amount of matter affects the matter term on the r.h.s. of the equation of 
motion ([10]), but the amount affected cancels out with the term (1 + a) in the denominator of (p!Q|) stemming from the 
new Newtonian potential. On the other hand, the factor Fsf augments (diminishes) for positive (negative) values of 
a for small distances compared to A, resulting in more (less) structure formation for positive (negative) values of a 
compared to the FCDM model. For r ^ X the dynamics is essentially Newtonian. 
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III. RESULTS 



A. Cosmological simulations of the Santa Barbara cluster 



In this section, we present results of the cosmological simulations of a FCDM universe with and without SF 
contribution. The initial condition of the system corresponds to the well known Santa Barbara cluster data that we 
get from the Heitmann's Cosmic Data Bank web page (http : //tSweb . lanl . gov/peo ple/heitmann/test3 . html[ ) . 
The initial condition uses a box 64 Mpc size, 128^ particles and start at redshift z = 63. 

In 1999 Frenk et al. [l| reported results of an extensive code comparison project involving twelve different codes. 
The aim of the project was to compare different techninques for simulating the formation of a cluster of galaxies 
(by now widely known as the Santa Barbara cluster) in a flat cold dark matter universe and to decide if the results 
from different codes were consistent and reproducible. We have repeated this test but restricted ourselves only to 
the dark matter component of the test. The cosmological parameters used in the simulation are as follows. We used 
128^ particles with initial positions and velocities with Qdm = 1, = 0, Qb = 0, Hq = 50 km/s/Mpc, box size = 
32 Mpc//i, (78 = 0.9 for the present-day linear rms mass fluctuation in spherical top spheres of radius 16 Mpc. The 
FCDM universe contents is given in units of the critical density, pc- The initial fluctuation spectrum was taken to 
have an asymptotic spectral index, n = 1, and the shape parameter F = 0.25, the value suggested by observations of 
large-scale structure [27]. See [1] and Heitmann's web page above for full details. In the shown snaps at z = 0, Fig. 
1, the big cluster near the center of the frame is the SB cluster. 

The particle masses are ^ 0.434 x 10-^^ Mq/H. The individual softening length was 50 kpc/h. These choices of 
softening length are consistent with the mass resolution set by the number of particles. 

We now present results for the FCDM model previously described. Because the visible component is the smaller 
one and given our interest to test the consequences of including a SF contribution to the evolution equations, our 
model excludes gas particles, but all its mass has been added to the dark matter. We restrict the values of a to the 
interval (—1, 1) [23, 24, 26] and use A = 1 Mpc//i, since this scale turns out to be an intermediate scale between 
the size of the clump groups and the separation of the formed groups. 

In Fig. 1 we show x-y snapshots at redshift 2; = of our FCDM model. Fig. 1 (a) presents the standard case 
without SF, i.e., the interaction between bodies is through the standard Newtonian potential. In (b) we show the 
case with a = 1, A = 1 Mpc/h. In (c) a = — 1/2, A = 1 Mpc/h. In (d) a = —1/4, A = 1 Mpc/h. One notes clearly 
how the SF modifies the matter structure of the system. The most dramatic cases are (b) and (c) where we have 
used a = 1 and a = —1/2, respectively. Given the argument at the end of last section, in the case of (b), for r <C A, 
the effective gravitational pull has been augmented by a factor of 2, in contrast to case (c) where it has diminished 
by a factor of 1/2; in model (d) the pull diminishes only by a factor of 3/4. That is why one observes for r < A more 
structure formation in (b), less in (d), and lesser in model (c). The effect is then, for a growing positive a, to speed 
up the growth of perturbations, then of halos and then of clusters, whereas negative a values (a ^ — 1) tend to slow 
down the growth. 

Next, we find the groups in the system using a friend-of- friend (FOF) algorithm and select one of the most massive 
ones. The chosen group is located approximately at the center in Fig. 1, the SB cluster. The group was analyzed by 
obtaining their density profiles (Fig. 2(a)) and circular velocities (Fig. 2(b)). The more cuspy case is for a = 1 and 
the less cuspy is for a = —1/2. The circular velocity curves where computed using v'^ = GNM{r)/r. The case with 
a = 1 corresponds to higher values of Vc, since this depends on how much accumulated mass there is at a distance r 
and this is enhanced by the factor Fsf for positive values of a. 



We finish this work by discussing a possible way to find a dark matter halo equation of state (EOS). The EOS may 
be useful to characterize more completely the state of a dark matter halo and could be a way to discriminate dark 
matter models. 

In general relativity both density and pressure contribute to modify the space-time geometry. For a static and spher- 
ically symmetric system in equilibrium like a dark matter halo we can not a priori neglect the pressure contribution. 
Therefore, and in the Newtonian limit we must solve the equation (see Misner et al. [2 81] ) 



B. Equation of state of a dark matter halo 




(12) 



together with the equation of hydrostatic equilibrium condition. 



dp . ^d^ 



(13) 
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FIG. 1: x-y snapshots at z = of a FCDM universe. See text for details. 




FIG. 2: (a) Density profiles of the SB cluster at z = 0. The cluster is located at the center in Fig 1. Vertical scale is in units 
of pQ — lO^°M0/i~^/(/i~^kpc)^. (b) The corresponding circular velocity. 
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EOS : NFW and Pseudo-isothermal : p-p'* 




10-5 0.001 0.1 10 

P 

FIG. 3: Equation of state for two density profiles. Pseudo-isothermal profile which has a power law behavior, p ^ with 
n = 1 and the NFW profile which has two power law behavior, one with n — lA and another with n — 1. Pressure and density 
are in geometrical units 



where p{r) = Pr{r) + 2pt{r) and we have assumed spherical symmetry. We also use the flatness condition on the 
circular velocity at large distances, 



2 _ ^ ^^N _ f^Q^^^^^^ 
dr 



to construct the boundary condition 



dp 
dr 



= -{p{Rh)^p{Rh)) 



r=RH 



Rh 



where Rh is the size of the system. The other condition is 



p{Rh) = 



(14) 



(15) 



(16) 



by, 



Results for two density profiles are given in figure HITBI The dashed line is for the isothermal density profile given 



piso 



1 + {r/vcY 



(17) 



where pc is the core density and Vc is the scale length of the matter distribution given by the isothermal profile. 
Whereas solid line is the Navarro-Frenk- White (NFW) density profile given by 



PNFW 



(r/r,)(l + r/r,)2 



(18) 



where pes is the NFW typical density and Vg is the scale length of the matter distribution given by the NFW profile. 

We may observe that NFW EOS has two power laws, i.e., behavior p ~ p^, for low densities n = 1.4 and for high 
densities n = 1 that is the behavior of the isothermal EOS. We may compare this results with the EOS of the SB 
cluster, shown in Fig. IIII B[ Where we have assumed in order to do the calculations, that p = pcr^, here cr^ is the 
radial dispersion of velocities of the cluster. This numerical EOS is shown in figure IIIIBi We may notice the power 
law behavior, n = 1 which corresponds to an isothermal density profile or to the high density case for the NFW EOS 



IV. CONCLUSIONS 



The theoretical scheme we have used is compatible with local observations because we have defined the background 
field constant < >= G^{1 + a). A direct consequence of the approach is that the amount of matter (energy) has 
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FIG. 4: The SB cluster extracted from the numerical simulation of the CDM model using a FOF algorithm. 

EOS : Santa Barbara cluster : p-^p", n=l 
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FIG. 5: Equation of state for SB cluster (dots), po is the unit of pressure 6.77 x 10 Pascals. The solid line is for an equation 
of state p ^ with n = 1. 



to be increased for positive values of a and diminished for negative values of a with respect to the standard FCDM 
model in order to have a fiat cosmological model. Quantitatively, our model demands to have a) = 1 and this 

changes the amount of dark matter and energy of the model for a flat cosmological model, as assumed. The general 
gravitational effect is that the interaction including the SF changes by a factor Fsrir^ X) ^ 1 -\- a (l + J) for r < A 
in comparison with the Newtonian case. Thus, for a > the growth of structures speeds up in comparison with the 
Newtonian case. For the a < case the effect is to diminish the formation of structures. For r > A the dynamics is 
essentially Newtonian. 

Additionally we have found numerically and EOS for the SB cluster. However, we assume that p = pa'^, where p 
and (Jr were obtain from SB cluster particle data. We leave for a future paper the numerical computation of Pr and 
Pt from the pressure tensor of the SB cluster and more detailed analysis. We have compared the power-law behaviour 
{p ~ p^) of the SB cluster EOS with the corresponding behavior of an EOS obtained solving Eqs. (p!2|) and (p!3|) for 
two density profiles, isothermal profile (pTj) and NEW profile (p!8|) . We see that the SB cluster is like the isothermal 
profile for some range of lower densities. Even more NEW profile has two power-law, for low densities behaves as 
n = 1.4 and for high densities as n = 1 that corresponds two the isothermal profile. 
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oration. The simulations were performed in the UNAM HP cluster Kan-Balam. 
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